First-principles modeling of multferroic RM^Os 
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We investigate the phase diagrams of RMr^Os via a first-principles effective-Hamiltonian method. 
We are able to reproduce the most important features of the complicated magnetic and ferroelectric 
phase transitions. The calculated polarization as a function of temperature agrees very well with 
experiments. The dielectric-constant step at the commensurate-to -incommensurate magnetic phase 
transition is well reproduced. The microscopic mechanisms for the phase transitions are discussed. 



RMn 2 5 (R=Tb, Dy, Ho, Y etc.) belong to a 
very special class of multiferroics, because the ferro- 
electricity is driven by the magnetic ordering [J, 0, 
These compounds therefore possess strong magnetoelec- 
tric (ME) coupling, showing remarkable new physical 
effects, such as the colossal magnetodielectric [4| and 
magneto-polarization- flop effects 0, 0, 0], etc. The 
strong ME coupling effects are not only interesting in 
the view of fundamental physics, but also they have po- 
tential important applications in future multifunctional 
devices. 

Because of the complex magnetic interactions and the 
ME coupling, RM^Os compounds undergo several mag- 
netic and associated electric phase transitions [il. [jl. r9L flG] 
upon cooling from room temperature to near zero tem- 
perature. Generally, these compounds transform at 
about 40 K from a paramagnetic (PM) phase to an an- 
tiferromagnetic (AFM) phase whose magnetic ordering 
is initially commensurate (CM) along the a axis. This 
phase transition is accompanied by a ferroelectriclike 
transition, with the appearance of a spontaneous po- 
larizations and a divergence of the dielectric constant. 
When the temperature is lowered further to about 20 K, 
the magnetic structures become incommensurate (ICM) 
along the a axis, and there is a drop of the electric po- 
larization together with the appearance of a step in the 
dielectric constant 0, 0> @] • The special phase transition 
sequence 0, flit is very puzzling and the driving forces for 
the phase transitions are not understood. It is therefore 
very important to explore the closely related magnetic 
and electric phase transitions to gain a full understand- 
ing of the microscopic mechanism of the ME coupling 
and novel physics in these materials. 

Recent neutron scattering experiments 11 , 12} as well 
as first-principles calculations [D, [l3| suggest that the 
strong ME coupling in RMn 2 Os is due to the "exchange 
striction" effect. However, previous first-principles calcu- 
lations [H, [H} were limited to zero temperature, and did 
not provide information about the phase transitions. The 
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FIG. 1: (Color online) A schematic sketch of the magnetic 
structure projected onto the ab plane. The angle of spin ori- 
entation is given in the parentheses. The dashed, solid, and 
double lines represent J3, J4 and J5 exchange interactions 
respectively. Ji and J2 (not shown) are along the c-direction. 



phase diagrams of RM^Os materials have been studied 
via a phenomenological approach This approach, 

based on symmetry considerations only, does not reveal 
any of the microscopic mechanisms of the ME coupling. 
In this letter, we present a first study of the phase dia- 
grams of RMn2 05 materials as a function of temperature 
by using a first-principles effective-Hamiltonian method 
[15j . We obtain the most important features of the phase 
diagram, including the the magnetic PM-CM-ICM tran- 
sitions, the accompanying ferroelectric transitions, the 
electric polarization as function of the temperature, and 
the dielectric-constant step at the CM-ICM transition. 

The high-temperature crystal structure of TbM^Os is 
orthorhombic (space group Pbam) with four TbM^Os 
formula units per primitive cell, containing Mn + 06 oc- 
tahcdra and Mn 3+ Os pyrami ds [l 61) - The effective Hamil- 
tonian was derived in Ref. 13 from a Heisenberg-like 
model. The spin-phonon coupling comes from the de- 
pendence of the exchange interactions J a on the phonon 
modes u\. J a {{u\\) was expanded around the high- 
symmetry structure to second order in the phonon mode 
amplitudes. Five nearest-neighbor (NN) exchange inter- 
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actions were included, as sketched in Fig. [TJ J3 is the 
Mn 4+ -Mn 3+ superexchange interaction through pyrami- 
dal base corners, while J4 is the superexchange interac- 
tion through the pyramidal apex ll| . The Mn 3+ ions in 
connected pyramids couple to each other antiferromag- 
netically through J5, whereas J\ and J2 couple Mn 4+ 
ions along the c axis. 

At lower temperature a further distortion occurs, re- 
ducing the crystal symmetry to Pb2im. The lattice dis- 
tortion involves 14 IR-active Bi u modes [H, U3]. Since 
the symmetry-lowering displacement is extremely small, 
we treated this displacement (henceforth u) as the only 
phonon normal mode in the model. Only the single pa- 
rameter J^=dJs/du was assumed to be involved in the 
first-order spin-phonon interaction. The neglect of 
terms, which renormalize the phonon frequencies and 
lead to the phonon anomalies near the magnetic phase 
transitions 17, 18, HH, is justified because these terms 
have a quite small effect on the phase diagrams studied 
here. The simplified Hamiltonian is then 



E ({ u k}) = E + ^2-muj 2 ul + ^2^ k iu k U[ (1) 



ije.J a ij£-h k 

Here Eq is the energy of the high symmetry structure 
without magnetic interactions, whereas m and lo are the 
reduced mass and frequency of the IR-active mode, u k is 
the A:-th local phonon mode , and are force-constant 
matrix elements that couple the NN local phonon modes. 
This last term was absent from Ref. [HI, but is included 
here to describe the phonon dispersion properly. We as- 
sume that the are isotropic in the ab plane, and we 
neglect the much smaller couplings along the c direction. 

Since the RM^Os compunds have similar phase di- 
agrams, we choose TbM^Os as an example, and de- 
termine the parameters for the simplified Hamiltonian 
Eq. (fT]) by carrying out a series of first-principles cal- 
culations on this compound [H, [ijj]. The calcula- 
tions were based on density-functional theory within 
the generalized-gradient approximation (GGA) imple- 
mented in the Vienna Ab-initio Simulations Package 
(VASP)[li Hj. Projector augmented- wave (PAW) 
pseudopotentials[22j and a 500 eV plane- wave cutoff were 
used. Spin polarization was included in the collinear ap- 
proximation. The resulting J parameters can neverthe- 
less be used to model noncollinear situations. 

To get the spin-phonon coupling constant J3, it is 
enough to use the the energy difference AE between 
the high-symmetry and the ground-state low-symmetry 
structures. To simplify the notation, we redefine u to 
be a dimensionless parameter taking the value of unity 
at the ferroelectric low-symmetry state, and assign spin 
moments |S,|=1.0 as well. Then it is easy to show that 
AE/4. We have J3 ~ 1.125 meV. In order to calcu- 
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late the phonon coupling constant we calculate the 
total energies of different local- mode configurations. In 
practice, we find that including the short-range phonon 
interaction only has a very small effect on the results. 
The exchange interactions J\ - J5 were fitted to the total 
energies of different spin configurations and were given 
in Ref. Il3l Alternatively, the exchange interactions can 
be calculated from the extended Kugel-Khomskii model 

In the present work, we have now also fitted 
the parameters to GGA+U calculations 24j] with 
1.0eV<C/<4.0eV on the Mn ions. We find that J' z de- 
creases with increasing U, falling to J' 3 — 0.325meV at 
U = 4.0 eV. As we shall see, this improves the com- 
parison of some of our later simulation results with ex- 
periment. Unfortunately, increasing U also worsens the 
agreement with experiment for the J parameters them- 
selves. This tension between the fitting of J and J' pa- 
rameters will be further discussed later. 

We investigated the finite-temperature behavior of our 
effective Hamiltonian by using Monte Carlo (MC) simula- 
tions. Traditional serial-temperature MC methods have 
great difficulty treating systems with complex frustrated 
interactions. Moreover, the present system has a first- 
order CM-AFM-to-ICM phase transition which would 
be very difficult to treat using conventional methods. 
Here we adopt the replica-exchange method (25| in which 
one simulates M replicas each at a different temperature 
T covering a range of interest, and allows configurational 
exchange between the replicas. Importantly, the inclu- 
sion of high-T configurations ensures that the lower-T 
systems can access a broad phase space and avoid be- 
coming trapped in local minima. 

We perform the simulations on an L x L x L cubic 
cell with periodic boundary conditions. Each unit cell 
contains eight spins and two local phonon modes. In 
the simulations, one MC sweep is defined to consist of 
a series of attempts of all variables. We performed the 
simulations at temperatures ranging from 3 to 90 K. The 
temperatures are adjusted to ensure that the exchange 
rates between adjacent remain close to 20%. At each T 
we carry out an initial 10 4 sweeps to prepare the sys- 
tem before allowing replica exchange. We discard these, 
as well as the first 10 6 sweeps after replica exchange is 
started, when computing equilibrium properties. Sample 
averages are accumulated over 2xl0 6 sweeps, without 
replica exchange to avoid a sign problem. 

We give here the results of typical simulations on a 
12x12x12 cell. Figures[2l[a) and (b) depict the polariza- 
tion and dielectric constant respectively, which are cal- 
culated via P—{u) and e—((u 2 ) — (u) 2 )/T. If we use 
J3 = 1.125 meV and the exchange interactions are fit- 
ted from the GGA calculations, we get a single magnetic 
PM-to-CM-AFM transition at about 58 K, accompanied 
by a ferroelectric transition (shown as the dotted lines in 
Fig. [J). This result misses the important CM-to-ICM 
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FIG. 2: (a) The electric polarization P as a function of the 
temperature; (b) The dielectric constant e as a function of 
temperature for both J3 = 1.125 and J3 = 0.4. The inset 
window shows the full view of e. The dielectric constants are 
normalized to unity at high temperature. 



phase transition and overestimates the PM-to-CM tran- 
sition temperature. The problem can be traced to the 
too-large spin-lattice coupling constant J3. Including the 
on-site Coulomb U can reduce J3, but at the same time 
it worsens the exchange interactions. We thus find that 
neither an effective Hamiltonian built on a pure GGA 
calculation, nor one built on GGA+U with a single value 
of U, can give good overall agreement with experiment. 

However, if we are willing to adjust the parameters by 
using the exchange coupling taken from pure GGA and 
the J3 from GGA+U, the situation improves dramati- 
cally. If J3 is reduced to about 0.4 meV, as obtained from 
GGA+U with U =3 eV, we obtain two phase transitions 
at about 42 and 18 K respectively, in very good agreement 
with experiment. The nature of each phase transition was 
identified via Fourier analysis of the spin configurations. 
At 80 K, the spins are fully disordered, indicating a PM 
phase. When T is lowered to 25 K, the spin spectrum 
shows a dominant peak at q = (0.5,0,0.5), suggesting a 
CM-AFM phase. The Fourier spectrum of the spin struc- 
ture at 5K shows dominant peaks at q=(5/12,0,0.5), in- 
dicating it is in the ICM phase. We therefore obtain the 
most important PM-CM and CM-ICM phase transitions, 
and the transition temperatures agree very well with the 
experimental values (~ 38-44 K for the AFM-CM order- 
ing along the a axis, and ^20 K for the ICM ordering) 
0, H S Eot The calculated q x =5/12 is slightly smaller 
than the experimental values (~0. 46-0. 48). It is worth 
noting that the calculated q x is restricted by the super- 
cell sizes in the simulation, which can be improved by 



increasing the supercell size. Simulations on a 14 x 14 x 14 
cell give q x =3/7. While we did not reproduce the correct 
q z in both the CM and ICM phases (probably because 
we only have NN interactions in the model Hamiltonian) , 
the fact that we nevertheless reproduce the correct phase 
transition sequence tends to confirm that the q z value is 
not important for the ME coupling in those materials 

The solid curves in Figs. EJa) and (b) show the spon- 
taneous polarization P and the dielectric constant e as 
functions of T for Jg=0.4 meV. The polarization increases 
strongly as the temperature is reduced through the PM- 
CM transition, but then it drops suddenly almost to zero 
at the CM-ICM transition. The magnetically induced 
polarization behaves as P <x (S3 ■ S4), where S3 and S4 
are the spins of the Mn 3+ and Mn 4+ ions coupled via 
the J3 interaction. In the ICM phase, S3 and S4 are 
almost orthogonal (i.e., 8 ~ 7r/2 in Fig. Q]), whereas in 
the CM phase they are parallel or antiparallel. These 
results are in excellent agreement with the experimental 
results for RMn 2 C)5 compounds 0, @, [ll|, [26[ , especially 
for YMn 2 5 [HI and HoMn 2 5 0. Note, however, that 
our simulation does not reproduce the reemergence of a 
polarized state observed experimentally in TbMn 2 Os at 
still lower temperature This is probably because we 
ignore the spins of Tb 4/ electrons in our model. Exper- 
imentally, it is observed that Tb is magnetically ordered 
below ~10K, which might play an important role in the 
reemergence of the polarization at low T [7j . 

The dielectric constant shows a peak at the PM-CM 
transition as a consequence of the ferroelectric phase 
transition. Most interestingly, the dielectric constant 
step at the CM-ICM transition has been well reproduced 
in the simulation, in which e jumps by about 75% in go- 
ing from the CM phase at 18 K to the ICM phase. The 
step is very interesting and important, because it may 
directly relate to the colossal magnetodielectric effect, 
which happens just at the CM-ICM transition tempera- 
tures in these materials. 

To gain a better understanding of the spin-lattice cou- 
pling effects on the magnetic and structural phase tran- 
sitions, we plot the J^-temperature phase diagram in 
Fig. [3] As we see, if the spin-phonon coupling is too 
strong (J3 > 0.42 meV), there is only a PM-CM transi- 
tion, and no CM-ICM transition (as in BiMn 2 5 [^V 
In contrast, if J3 is very small (J3 < 0.30 meV), the CM 
state will not appear, and instead a state having spin- 
glass (SG) character will appear above the ICM state. 
The complex nature of the phase diagram is due to the 
frustration of the J3 interactions. It is easy to see that 
the CM states do not have the lowest magnetic energies, 
since J3 induces spins to rotate to decrease the energy. 
Since the spins have the same wave vector q z = 0.5 along 
the c axis in both the CM and ICM phases, the inter- 
actions due to Ji and J 2 do not change in the the two 
phases, and can be neglected in the discussion. Accord- 
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FIG. 3: The J^-temperature phase diagram of the system, 
where PM, CM, ICM and SG represent paramagnetic, com- 
mensurate, incommensurate and spin-glass-like phase respec- 
tively. 



ing to the phase factors shown in Fig. [1] the energy of 
the ICM state can be written as 

Eicm ~ 8J4 — 2 J 5 cos(2nq x ) + 4 J 3 [cos 9 + cos(27rg x — 0)] . 

(2) 

Here we assume that two spins connected by J4 are al- 
ways antiparallel to each other, because the J4 interac- 
tions, each having two NNs, are much stronger than the 
J3 and J5 interactions. We also ignore the phonon con- 
tribution, because in the ICM phase, 9 ~ ir/2 and the 
energy from spin-lattice coupling J3 is small. For the 
CM state, we have 9 = 0, q x — 0.5. Therefore, 



E, 



CM 



8 J 4 + 2 J 5 - 4 J', 



(3) 



The energy difference between the CM and ICM phases 
is determined by the competitions among J3, J5 and J3. 
In the case of J3=1.125meV, the energy of the CM state 
is always lower than that of the ICM state, and there is 
no CM-ICM transition. However, when J3 decreases to 
0.4 meV, one can find suitable 9 and q x that allow the 
ICM state be the ground state. Since there is no group- 
subgroup symmetry relation between the CM and ICM 
states, the phase transition between them is necessarily 
a first-order one. We speculate that the transition occurs 
because the entropy of the CM state is larger than that 
of the ICM state for suitable J3. 

The above results are calculated from the L=12 cell. 
We have also obtained similar results for the L=10 and 
14 cells. However, due to the subtle nature of the ICM 
state, the parameters J3 and J3 have to be slightly ad- 
justed to produce results that are in good agreement with 
experiments for different cell sizes. 

The "semi-empirical" philosophy we have adopted here 
has been to start with first-principles derived parameters, 
make the minimal empirical modifications to the param- 
eters to get good agreement with experiment, and then 
use the resulting model to make predictions. The fact 



that we have to adjust J3 by hand (through the choice of 
U) to obtain good agreement with experiment is clearly 
somewhat unsatisfactory. However, there is considerable 
precedent for such an approach. For example, in simula- 
tions of ferroelectrics via effective-Hamiltonian methods, 
it is a common practice to adjust the lattice constant 
to agree with experiment through the application of a 
fictitious negative pressure (l5j . 

In summary, we have investigated the phase diagrams 
of RMn2 05 using a first-principles effective Hamilto- 
nian method. We obtained the most important fea- 
tures of the phase diagrams of multiferroic RM112O5 
compounds, including the sequence of magnetic and fer- 
roelectric phase transitions. Most importantly, we ob- 
tained the dielectric-constant step at the commensurate- 
to-incommensurate magnetic phase transition, which is 
key to understand the colossal magnetodielectric effects. 
The work further clarified the microscopic mechanism of 
the magnetoelectric coupling in RM112O5, and can be use- 
ful for exploring other multiferroic materials. 

L.H. acknowledges the support from "Hundreds of Tal- 
ents" program from CAS, and NNSF of China, Grant 
No. 10674124. D.V. acknowledges the support from NSF 
Grant DMR-0549198. 



[1 
[2 

[3. 
[4 
[5 
[6 

ir 

[8 
[9 
[10 

[11 
[12. 

[13 
[14 
[15 
[16 

[17 

[is; 

[19 
[20 
[21 
[22 
[23 
[24 

[25 
[26 
[27 



W. Eerenstein et al, Nature 442, 759 (2006). 
S.-W. Cheong et al, Nature Materials 6, 13 (2007). 

C. Wang et al, Phys. Rev. Lett. 99, 177202 (2007). 
N. Hur et al, Phys. Rev. Lett. 93, 107207 (2004). 
T. Kimura et al, Nature 426, 55 (2003). 

T. Goto et al, Phys. Rev. Lett. 92, 257201 (2004). 

N. Hur et al, Nature 429, 392 (2004). 

L. C. Chapon et al, Phys. Rev. Lett. 96, 097601 (2006). 

D. Higashiyama et al, Phys. Rev. B 72, 064421 (2005). 
H. Kimura et al, J. Magn. Magn. Mater. 321, 854 (2009). 
L. C. Chapon et al, Phys. Rev. Lett. 93, 177402 (2004). 
P. G. Radaelli et al, Phys. Rev. B 79, 020404(R) (2009). 
C. Wang et al, Phys. Rev. B 77, 134113 (2008). 

A. Harris et al, Phys. Rev. Lett. 100, 217202(R) (2008). 
W. Zhong et al, Phys. Rev. Lett. 73, 1861 (1994). 
J. A. Alonso et al, J. Phys.: Condens. Matter 9, 8515 
(1997). 

A. F. Garciia-Flores et al, Phys. Rev. B 73, 104411 
(2006). 

J. Cao et al, Phys. Rev. Lett. 100, 177205 (2008). 
T. Shen et al, Phys. Rev. B 78, 134413 (2008). 
G. Kresse et al, Phys. Rev. B 47, RC558 (1993). 

G. Kresse et al, Phys. Rev. B 54, 11169 (1996). 
P. E. Blochl, Phys. Rev. B 50, 17953 (1994). 

H. Das et al, Phys. Rev. Lett. 100, 186402 (2008). 

G. Giovannetti et al, Phys. Rev. Lett. 100, 227603 
(2008). 

R. H. Swendsen et al, Phys. Rev. Lett. 57, 2607 (1986). 

I. Kagomiya et al, Ferroelectrics 286, 167 (2003). 
A. Munoz et al, Phys. Rev. B 65, 144423 (2002). 



